source("/users/brucedesmarais/Dropbox/SenPE/Code/CreateSenPENets.R")

## Add centrality to Fowler's Data ##
id_master[,5] <- procName(id_master[,5])
fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")
fowid <- paste(fowler$labels,fowler$congress,sep="")

nc <- rep(0,nrow(fowler))
ec <- rep(0,nrow(fowler))
bc <- rep(0,nrow(fowler))
dc <- rep(0,nrow(fowler))

ncl <- rep(0,nrow(fowler))
ecl <- rep(0,nrow(fowler))
bcl <- rep(0,nrow(fowler))
dcl <- rep(0,nrow(fowler))

l_betweeness <- rep(NA,nrow(fowler))
l_connectedness <- rep(NA,nrow(fowler))
l_evcent <- rep(NA,nrow(fowler))
l_closeness <- rep(NA,nrow(fowler))
for(i in 1:length(wnets)){
library(igraph)
net <- wnets[[i]]
gnet <- graph.adjacency(net,mode="undirected")
congi <- as.numeric(names(wnets)[i])
distnet <- graph.adjacency(1/net,mode="undirected",weighted=T)
sp <- shortest.paths(distnet,algorithm="dijkstra")
csp <- c(sp)
csp_inf <- which(csp==Inf)
csp[csp_inf] <- 1.25*max(csp[-csp_inf])
msp <- matrix(csp,nrow(sp),nrow(sp)) 
new_cent <- 1/apply(msp,2,sum)
ev_cent <- evcent(gnet)
bet_cent <- betweenness(gnet)
cl_cent <- closeness(gnet)
pe_now <- colnames(net)
sm_now <- id_master[match(paste(pe_now,congi,sep=""),paste(id_master[,5],id_master[,1],sep="")),3]
sum_now <- id_master[match(paste(pe_now,congi,sep=""),paste(id_master[,5],id_master[,1],sep="")),2]
fow_now <- match(paste(sum_now,congi,sep=""),paste(fowler$labels,fowler$congress,sep=""))

nc[fow_now[which(!is.na(fow_now))]] <- new_cent[which(!is.na(fow_now))]
ec[fow_now[which(!is.na(fow_now))]]<- ev_cent$vector[which(!is.na(fow_now))]
bc[fow_now[which(!is.na(fow_now))]] <- bet_cent[which(!is.na(fow_now))]
dc[fow_now[which(!is.na(fow_now))]] <- (apply(net>0,2,sum)/2)[which(!is.na(fow_now))]

fow_now[which(is.na(fow_now))] <- 1
betweeness <- fowler$betweeness[fow_now]
connectedness <- fowler$connectedness[fow_now]
evcent <- fowler$evcent[fow_now]
closeness <- fowler$closeness[fow_now]


if(i <length(wnets)){
	sum_next <- id_master[match(paste(sm_now,congi+1,sep=""),paste(id_master[,4],id_master[,1],sep="")),2]
	fow_next <- match(paste(sum_next,congi+1,sep=""),paste(fowler$labels,fowler$congress,sep=""))
	ncl[fow_next[which(!is.na(fow_next))]] <- new_cent[which(!is.na(fow_next))]
	ecl[fow_next[which(!is.na(fow_next))]]<- ev_cent$vector[which(!is.na(fow_next))]
	bcl[fow_next[which(!is.na(fow_next))]] <- bet_cent[which(!is.na(fow_next))]
	dcl[fow_next[which(!is.na(fow_next))]] <- (apply(net>0,2,sum)/2)[which(!is.na(fow_next))]
	
	l_betweeness[fow_next[which(!is.na(fow_next))]] <- betweeness[which(!is.na(fow_next))]
	l_connectedness[fow_next[which(!is.na(fow_next))]] <- connectedness[which(!is.na(fow_next))]
	l_evcent[fow_next[which(!is.na(fow_next))]] <- evcent[which(!is.na(fow_next))]
	l_closeness[fow_next[which(!is.na(fow_next))]] <- closeness[which(!is.na(fow_next))]
	}

}

fowler <- data.frame(fowler,nc,ec,bc,dc,ncl,ecl,bcl,dcl,l_betweeness,l_connectedness,l_evcent,l_closeness)
fowler <- subset(fowler,is.element(fowler$congress,97:105))

# Scatter plot of centrality scores
for(i in 97:105){
filei <- paste("/users/brucedesmarais/Dropbox/SenPE/Tex/twocents",i,".pdf",sep="")
nci <- subset(fowler$nc,fowler$congress==i)
mnci <- median(nci,na.rm=T)
connecti <- subset(fowler$connectedness,fowler$congress==i)
mconn <- median(connecti,na.rm=T)
cent_rho = cor(nci,connecti,use="complete.obs")
pdf(filei,height=4.5, width=5, family="Times", pointsize=13)
par(las=1,mar=c(6,6.25,.5,.6),cex.lab=2,cex.axis=1.75)
plot(connecti,nci,pch=4,col="grey35",xlab="",ylab="",cex=1.35)
abline(h=mnci,v=mconn,lwd=2)
title(xlab="Cosponsorship",line=3.5)
title(ylab="Press Events",line=4.75)
rr <- round(cent_rho,digits=3)
legend("topleft",legend=substitute(paste(rho," = ",rr),list(rr=rr)),col="white",cex=1.5,box.col="white",inset=-.05)
dev.off()
}

# Simple Means and Tests of Hypotheses
library(sciplot)
med_nc <- median(fowler$nc)
med_ncl <- median(fowler$ncl)
ncid <- as.factor(fowler$nc > med_nc)
nclid <- as.factor(fowler$ncl > med_ncl)

filei <- "/users/brucedesmarais/Dropbox/SenPE/Tex/barplot.pdf"
pdf(filei,height=4, width=5, family="Times", pointsize=13)
par(las=1,mar=c(4,4,2,.6),xpd=T)
bargraph.CI(ncid,fowler$pb,nclid,xaxt="n",legend=T,leg.lab=c(expression(paste("Low Centrality ", italic(t))),expression(paste("High Centrality ", italic(t)))),xpd=T,x.leg=1,y.leg=10.25,ylab="Mean Bills Passed",col=c("grey45","grey80"))
axis(1,at=c(2,5),labels=c(expression(paste("Low Cent. ", italic(t-1))),expression(paste("High Cent. ", italic(t-1)))))
dev.off()

write.dta(fowler,"fowler_plus.dta")
